Multiple phase slips phenomena in mesoscopic superconducting rings 
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We investigate the behavior of a mesoscopic one-dimensional ring in an external magnetic field by 
simulating the time dependent Ginzburg-Landau equations with periodic boundary conditions. We 
analyze the stability and the different possible evolutions for the phase slip phenomena starting from 
a metastable state. We find a stability condition relating the winding number of the initial solution 
and the number of flux quanta penetrating the ring. The analysis of multiple phase slips solutions 
is based on analytical results and simulations. The role of the ratio of two characteristic times u is 
studied for the case of a multiple phase slips transition. We found out that if u >> 1, consecutive 
multiple phase slips will be more favorable than simultaneous ones. If u « 1 the opposite is true 
and we confirm that u >> 1 is often a necessary condition to reach the ground state. The influence 
of the Langevin noise on the kinetics of the phase transition is discussed. 



I. INTRODUCTION 

Nonequilibrium phenomena in superconductors are a 
challenging area for both experimental and theoretical 
research. They are also crucial for the development 
of applications as they are the key to the appearance 
of resistive states in superconducting samples and to 
the possible use of vortex dynamics. In one-dimension, 
the resistive phase slip process foreseen by Little^ has 
been described quantitatively by the theory developed by 
Langer and Ambegaokar— and extended by McCumber 3 
and Halperini (LAMH). The LAMH theory describes 
thermally activated phase slips. It evaluates the resis- 
tance of a ID superconductor when driven out of thcr- 
modynamical equilibrium by a voltage or current source. 
Since then, this theory has been accepted in rather good 
agreement with experiments (for a more complete review 
of theoretical and experimental works, see Ref. [5|). 

More recently, simulations were carried out on different 
versions of the time-dependent Ginzburg-Landau (GL) 
(TDGL) equations in order to investigate the dynamics 
of the process. In particular, the case of the ID ring 
has raised interest since it exhibits multiple metastable 
states which can be reached by phase slip processes. In 
the work of Tarlie and Elder—, the ring was submitted to 
an electromotive force which constantly accelerates the 
superconducting electrons. Another approach was used 
by Vodolazov and Peeters^, who increased the magnetic 
field gradually. The latter simulations were confirmed by 
the experiments made in Refs. H and[^. 

In the present work, we consider a superconducting 
ring of thickness d, radius R and length L as represented 
in Fig. Q] For simplicity, we consider d « £ << A c ff, 
R > £ and R « A e ff, where £ is the coherence length 
and A e ff is the PearliS penetration depth. The first two 
conditions allow us to treat the ring as one dimensional 
and the last two conditions account for the mesoscopic 
size of the ring. We apply an external constant magnetic 
field H, perpendicular to the ID ring. The field H is an 
external parameter that controls the parameters (super- 
fluid density, current etc.) of the ring. 

Our approach is different from that of Refs. @ and0 be- 



FIG. 1: Geometry of the ring of radius R and thickness d with 
the magnetic field H. 



cause we investigate the stability of a stationary state and 
observe the relaxation process from the initial metastable 
state to a stable one, without any evolving external pa- 
rameter influencing the dynamics. Indeed, starting from 
a solution of the stationary GL equations, we derived the 
stability condition regarding small time-dependent per- 
turbations. For example starting from a stable state, an 
increase in the magnetic field can result in the supercon- 
ducting state to evolve toward a new equilibrium. This 
new equilibrium is found via one or more phase slip pro- 
cesses. Moreover, using an instant increase in the mag- 
netic field in contrast with the ramp applied in Ref. l7| 
allows us to investigate the competition between consec- 
utive and simultaneous phase slips. 



II. THE TIME DEPENDENT EQUATIONS 

The TDGL equations have been derived in different 
ways in order to describe the different nonequilibrium 
properties of superconductors. The range of validity of 
the TDGL equations has been widely discussed in the 
literature (see, for example, Refs. fTTI andfl^). It requires 
not only slow variation in the order parameter in space in 
comparison with the diffusion length but also slow vari- 
ation in the order parameter in time in comparison with 
the inelastic scattering time. Therefore the TDGL equa- 
tions provide an accurate description of kinetics if the size 
of the sample is larger or comparable with the coherence 
length £ and if the time scale is larger than the inelastic- 
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scattering time. Here we use the simplest version that 
allows to study the phase slip phenomena. It takes into 
account the presence of magnetic field and the possibility 
of charge imbalance by using the vector potential defined 
in the symmetric gauge and the electrostatic potential as 
introduced by Gor'kov and Kopnin.— The first equation 
in dimensionless units reads as: 



*-*|*| 2 -(iV + A) 2 * + ?7. (1) 



The dimensionless order parameter ^ = p(r, £)e ie ( r,t - ) 
takes the value \& = 1 in the equilibrium and in the ab- 
sence of field. The space and time coordinates (r, t) are 
measured in units of the coherence length £ and of the 
characteristic relaxation time of the phase Tg = 47rA rff°" 
(<r„ is the conductivity of the normal state and c is the 
speed of light). The Pearl penetration depth A c ff = 4- 
has been used instead of the London penetration depth A 
since the thickness d is small. The vector potential A is 
written in units of (4> Q is the flux quantum) and the 

electrostatic potential <& in units of , with e being the 
electron mass and h as the reduced Planck constant.— 

The only dimensionless parameter left in the equation 
is the ratio u — — between the two characteristic times: 

re 

t p is the characteristic time of the evolution of the ampli- 
tude of the order parameter, whereas Tg accounts for the 
dynamics of the phase. The estimates from microscopic 
theories give the value of u ranging from 5.79 to 12 (see 
Refs. GJandQJ). However, we consider here the general 
case where < u < oo. 

We also add a Langevin noise rj during the simulations. 
The intensity of the noise may vary considerably from 
one material to another and depends on the experimental 
conditions. Nevertheless, for the case of the second-order 
phase transition, the noise should be small and on the 
order of Jr§r-z ~ y-(%^) 2 << 1 in dimensionless units 

when d < £. Here, ks is the Boltzman constant, T is the 
temperature, T c is the critical temperature, H C2 is the 
second critical field and Ef is the Fermi energy. 

The first TDGL Eq. (fl]) is very similar to what is some- 
times referred to as the complex GL (CGL) equation, 
except for the Langevin noise and the terms containing 
electrostatic potential and vector potential. The second 
TDGL equation is the decomposition of the total current 
into the superconducting and the normal currents 



V x (V x A) = j 
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The total current j = V x (V x A) is the sum of the su- 
perconducting current j s = -|(\E'*V , I' - WW*) - A|*| 2 
and the normal current j n = — 4^ — V$ in dimension- 
less units. The GL parameter k = is large and the 



dimensions of the ring are small. Therefore, we neglect 
the corrections to the vector potential. Moreover, using 
the electroneutrality relation divj = (see Ref. LLl| for 
details), Eq. is simplified as 



-(W*V* - WW*) + A|*| 



(3) 



The stability of the stationary solutions of the CGL 
equation has been studied in details (see Ref. [l6h . and 
in ID, one can derive the conditions of Eckhaus-type in- 
stability. In the presence of magnetic field, most of the 
results are comparable. Indeed, we can rewrite Eq. ([1]) 
in the same way as in Ref. |l_7j by separating real and 
imaginary parts and neglecting the noise: 
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(5) 



where p is the amplitude, 9 is the phase of the com- 
plex order parameter 5', and x is the longitudinal dimen- 
sionless spatial coordinate. The superconducting current 
/0 2 (f§ — A) needs to be uniform for stationary solutions. 
It corresponds to zero electrostatic potential and leads to 

We use a similar derivation of the constant solutions 
as in Refs. [H and [l?] by describing Eqs. ((4]) and ([5]) as 
the equations of motion of a particle in a potential. We 
obtain a family of solutions in the well-known twisted 
plane-wave form, 



* fe = y/1 - (A - kf 



J(kx+9 ) 



(6) 



where k and 9q are real numbers. Periodic boundary con- 
ditions imply that k — ^j^-, where I is the dimensionless 
length of the wire and n is an integer that corresponds to 
the winding number or vorticity of the solution. Indeed, 
we have n — § || . These solutions are only possible 
when j s < j c — 2/\/27, j c being the GL critical current. 

The stability of such solutions has already been stud- 
ied in different cases, mostly in the context of the CGL 
equation (without including electromagnetic fields). In 
our case, we include in the analysis both the vector po- 
tential and the electrostatic potential. We then ana- 
lyze the stability of a solution in the form of Eq. ^ 
by linearizing the equations disturbed by a perturbation 
y{x, t) = yr{x, t)+iyi{x, t) and performing Fourier trans- 
form y r (q, t) + iyi(q, t). We obtain the system of algebraic 
equations, 
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y\\ _ i ( -2[l-{k-A) 2 )]-q 2 2i(k-A)q \/y\ 

y\J ~" {-2 l (k-A) q ~ 2 -^^[i-(k^A) 2 ] -q 2 -u[i-{k-A) 2 ])\ Vl )■ W 

The eigenvalues are consistent with the ones found with neither electrostatic potential nor magnetic field 6 and in the 
two-dimensional case J£ Only one of the two families of eigenvalues of this system can be positive, 

A = - y/A - Au + u 2 + (-8 + IQq 2 + 24u - 2u 2 ){k - A) 2 + (4 - 20u + u 2 ){k - A) 4 

-l-q 2 - l -u+ l -u{k-Af + {k-Af. (8) 

I 



Therefore, the stability condition reads as 



(9) 



For the case of periodic boundary conditions such as in 
the ID ring, we use the Fourier expansion of the pertur- 
bation, 



y(o) = c ™(y) s [q- y) » 



(10) 



where c m (y), m € Z are the Fourier coefficients of y and 
S is the delta function. Therefore, the stability condition 
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is in agreement with those found in previous works. 6,7 
Let us rewrite this condition in the form, 



FIG. 2: Shape of the free energy of the ring depending on the 
ratio and the winding number n. 



> l 2 [6(k - A) 2 ~ 2] . 



(12) 



It means that the higher modes corresponding to multiple 
simultaneous phase slips will be stable if the length of 
the wire is small enough. In other words, the length of 
the wire can restrict the number of phase slips that can 
happen at the same time. Moreover, the size of the ring 
plays an important role in the selection between multiple 
and consecutive processes as we discuss in Sec. IIII CI 
Depending on the size of the ring, different modes can 
have highest eigenvalue. 

We simplify the stability condition in: 
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(13) 



This condition may be viewed as a generalization of Eq. 
(jTTJ). It corresponds to the thermodynamically stable 
supercurrent j s = p 2 \J\ — p 2 (see, for example, Fig. 18. 
in Ref. [l2h . Rewriting the condition (|13p in dimensional 
units, we find a stability condition relating the winding 
number n of the solution with the number 4- of flux 
quanta penetrating the ring: 



< 



R 



(14) 



This condition is consistent with the ground state found 
by minimizing the free energy with a solution in the form 
of Eq. ©. As seen in Fig. O the ground state is reached 



when the value of 



is minimal. 



III. 



PHASE SLIP SIMULATIONS 



A. Mathematical formulation of the problem 



According to Eq. (|13p. the winding number needs to 
be changed in order to reach a stable state. The new so- 
lution will be closer to the ground state. Therefore, the 
superconductor should reach a state with lower current 
and energy. Indeed, remembering that the superconduct- 
ing current is j s = -A|*| 2 - i(**V* - W**), the 
transition to a state of lower current cannot occur with- 
out modifying the phase and hence the value of k. We 
observe a transition from the solution to a different 
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solution ^fej. Because of periodic boundary conditions, 
the transition from k to k\ requires discontinuity in the 
phase. Since the order parameter is a single- valued func- 
tion, this discontinuity or jump should be of 2irn, n being 
an integer number. Moreover, 'J = at the center of the 
phase slip to fulfill the continuity condition for the order 
parameter. Therefore, the transition between two solu- 
tions is a resistive phase slip. Once the slip is achieved, 
the phase recovers continuously. 

In our simulation, we start with k — 0. We expect 
the evolution to a state k\ = after a number rt\ of 
2tt phase slips (or an equivalent number of greater phase 
slips) corresponding to the final phase: 



27rni 



(15) 



If the number n\ of phase slips is equal to the integer 
part of the ratio 4^ , the final state reached is the ground 
state, but as we show, this is not always the case and it 
is difficult to predict (see Sec. MI CI for details). In other 
words, we observe the transition after an infinitesimal 
perturbation from the solution 



to another solution after n\ phase slips 



(16) 



^ 2.^ = y/l- (A-2irn 1 /l) 2 e^ 27lnix/l+e ° ) . (17) 

This solution may correspond to the ground state or a 
new metastable state of lower energy. 



B. The single phase slip 

In the simulations, we use fast Fourier transforms for 
spatial derivatives and Runge-Kutta method for time 
evolution. As discussed in Sec. IIII Al and as described 
in the LAMH theory, the amplitude of the order parame- 
ter vanishes in a narrow region near the phase slip center 
for a short time (see Fig. 19 in Ref. Il2h . Then the ampli- 
tude relaxes to a uniform state with a higher value than 
it had initially. The phase of the order parameter devel- 
ops a sharper sinusoidal form, until the minimum and the 
maximum disconnect, at the same moment when the am- 
plitude vanishes and in the same region, as discussed in 
Ref. [Ht Afterward, it relaxes to a sawtooth pattern cor- 
responding to the new state (fTS"]) so that both amplitude 
and phase are consistent with the solution [)17p. 

The superconducting current behaves similarly to the 
amplitude of the order parameter. It vanishes at the 
same point and at the same time as the phase slips. This 
region becomes resistive at that moment. Contrary to 
the amplitude, the process occurs with a decrease in the 
supercurrent, as explained in Sec. IIII Al 

According to the Josephson equation = A$, the 
electrostatic scalar potential and the phase are directly 



connected (A6 and A$ are respectively, the phase and 
electrostatic potential difference taken between two arbi- 
trary points). As the amplitude of the order parameter 
goes to zero, the average speed of the electrons is re- 
duced and thus some of them gather before the region 
of the phase slip. Therefore, there is a deficit of elec- 
trons at the end of the region and the electric field is 
created. This voltage influences the phase of the order 
parameter, until it makes a slip of 2ir when the amplitude 
vanishes. Then, the amplitude recovers and the voltage 
relaxes back to zero, as the electrons spread again uni- 
formly along the wire. The "discontinuity" in the phase 
remains (see Fig. 3 in Ref. ). 

The characteristic values for the process are deceiving 
because they depend strongly on the material and the 
thickness of the ring. For example in NbN we estimate 
Tg w lps. The total time of the transition between two 
states separated by a single phase is thus of the order 
of lOOps and the characteristic voltage appearing in the 
ring is of about 10//V. 



C. Multiple phase slips solutions 

Here we investigate the case of a transition where the 
stable state and the initial state are separated by more 
than one phase slip. The stability condition (| 13[) is insuf- 
ficient to predict the final state and the number of phase 
slips that occur. The best prediction of the number of 
phase slips is given by finding the winding number of the 
ground state (see Fig. [2]) . In general, the maximum num- 
ber of phase slips that can happen is equal to n — . 

Starting with n = 0, a high magnetic flux will thus be a 
good condition to observe multiple phase slips. However, 
we found out that the ground state is reached only for 
certain values of the parameter u. 

If u » 1, the phase 9 relaxes faster than the ampli- 
tude p. It favors multiple phase slips to happen at the 
same spot, one after another, as seen in Fig. |3fa). The 
oscillations of the amplitude of the order parameter dur- 
ing the time of the transition are similar to that observed 
in Ref. 0- In this case, the order parameter almost al- 
ways relaxes to the ground state, even if an intermediate 
state is stable according to condition (TT3"|) . 

If u « 1, the amplitude p relaxes faster than the 
phase. It favors processes where multiple phase slips 
happen more or less simultaneously at different phase 
slip centers as seen in Fig. EJb). Therefore, after a cer- 
tain number of simultaneous phase slips, the amplitude 
relaxes to a new uniform state which will be the final sta- 
ble state. In that case, the final state is not necessarily 
the ground state. These observations are consistent with 
the observations made in Ref. [?] that for u « 1 the final 
state may be "further away" from the ground state than 
for u » 1. 

Mathematically, the importance of the parameter u is 
emphasized by plotting different eigenvalues as a function 
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FIG. 3: Evolution of the distribution of the order parameter amplitude p for two consecutive phase slips (a) and two simultaneous 
phase slips (b). 



of the size of the ring. For u » 1, the highest eigenval- 
ues always corresponds to single phase slips happening 
one at a time. However, for u << 1, different eigen- 
modes have highest eigenvalue depending on the size of 
the ring. For large sizes, the simultaneous multiple phase 
slips modes are the most unstable (see Fig. 

It is possible to differentiate three cases when multi- 
ple phase slips happen: (1) there is only one phase slip 
center and the amplitude of the order parameter reaches 
zero several times at the same location before relaxing 
to a stable state. This is the general consecutive phase 
slips case. (2) There are more phase slip centers and the 
amplitude reaches zero at the same time at all locations. 
(3) There are more phase slip centers but the amplitude 
reaches zero at a different time for each location. 

Here we consider cases (2) and (3) as simultaneous 
phase slips. Indeed, the amplitude starts to decrease at 
the same time for all phase slip centers, so that the pro- 
cesses are called simultaneous. In Fig. OUa) one can also 
see the competition between simultaneous and consec- 
utive phase slips. Indeed, in that simulation, the am- 
plitude decreases first on two separate spots, but in the 
end, both phase slips happen consecutively on the same 
spot. Figure [5] shows the different behaviors depending 
on the parameter u. We confirm that the ground state 
is only reached for u >> f. Indeed, for those simula- 
tions, we have 4^ ~ 5 and we observe five phase slips 
only when u = 100. However, the number of phase slips 
does not necessarily grow monotonically with the value 
of u. Indeed, for u = 1, we observe less phase slips than 
for u = 0.5 and u = 0.1. 

According to Eq. (fT2"|) and Fig. [H the size of the 
ring restricts the number of phase slip centers. However, 
at small sizes, multiple phase slips are often ruled out 
because the magnetic flux is too low or too high. In 
particular, the magnetic flux induces more phase slips 
than the number of phase slip centers allowed by the 
eigenmodes if 



This condition restricts greatly the choice of parameters. 
In our case, with k = 0, we need A < 
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(18) 



-j=. Therefore, it limits the simulations to an area which 
is almost stable and therefore intrinsically not prone to 
simultaneous phase slips. 

As described in Eq. ([1]), we add a Langevin noise. The 
amplitude of the noise corresponding to typical experi- 
mental values is on the order 10 -3 which does not create 
any noticeable difference with the noiseless situation (all 
simulations excepting those shown in Fig. [5] were per- 
formed using a Langevin noise of this amplitude) . How- 
ever, according to our simulations the increase in the level 
of the noise leads to a stronger instability. We notice that 
the phase slips are accelerated, especially in the begin- 
ning of the process. In the multiple phase slips case, 
the behavior of the order parameter can be significantly 
altered. The phase slips tend to happen "less simultane- 
ously" and we even observe a mixing of consecutive and 
simultaneous phase slips [see Fig. (Hd)]. Rather surpris- 
ingly, the increase in the noise does not necessarily lead 
to the relaxation of the order parameter closer to the 
ground-state value. Sometimes we observe the opposite. 
For u = 0.5, increasing the noise reduced the number of 
phase slips from four to three as one can see by compar- 
ing Figs. EJb) and[|3b). We also note that the general 
conclusions of the stability analysis and the role of the 
parameter u stay valid. Last but not the least, in the 
absence of noise, a nonuniform perturbation needs to be 
applied in order to initiate the phase slip process. 

As pointed out in Ref. |7|, a strong magnetic field might 
trigger strong normal currents and destroy superconduc- 
tivity because of heating effects. In our simulations, a 
strong magnetic field is indeed required to trigger si- 
multaneous phase slips. However, in the u << 1 case, 
we observe simultaneous multiple phase slips already at 
magnetic fields on the order of 1G for an Al sample of 
length L — 40£ and £ = fOOnm. This magnetic field is 
smaller than the value described as the maximum field 
sustainable by such a ring in Ref. 8. The estimate of the 
heating of the ring may be made using Ohm's law, as we 
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FIG. 4: Eigenvalues A corresponding to different modes m as a function of the size I of the ring for u = 0.1 (a) and u — 100 
(b). Here k — and A — 0.8. When u » 1 there is no competition between the modes and consecutive phase slip processes 
will always dominate. With u << 1, the highest eigenvalue depends on the size of the ring. For larger length the simultaneous 
phase slips will dominate. 



consider that during a phase slip event, a part of the ring 
of length £ behaves like a normal conductor of conduc- 
tivity a n during a time Tg and sustains a current J s (in 
cgs units). The energy dissipated per unit of volume is 
then: 



E 



Jl [foA(l-A 2 )]' 
a n T6 16^ 2 A c 2 ff 



(19) 



This energy is less than a quarter of the condensation 

H 2 

energy 167 ^ 2 and such a small heating should not have 
any significant effect on the experimental detection of the 
phase slip. The case of multiple phase slips is more com- 
plicated and depends on the quality of the contacts of the 
sample with the thermostat. We believe that for the case 
of a free standing ring with a reasonable number of con- 
secutive phase slips, the local temperature does not reach 
T c because of the large time scales involved [see Figs. [Sfd) 
and[5je)]. In the case of simultaneous phase slips, simu- 
lations confirm that the phase slip centers spread along 
the whole sample and therefore the distances are large in 
comparison with the heat diffusion length [see Figs. E^a) 
and^b)]. In that case, the local heating is similar to the 
single phase slip case. 



IV. CONCLUSION 



We formulated the stability condition for the supercon- 
ducting state of a ID ring in a constant magnetic field. 
Using the TDGL equations, we found out that the state 
is stable when the difference between the vorticity and 
the number of flux quanta enclosed in the ring is small. 
The relaxation towards the stable state must therefore 
involve one or more phase slips. The study of the multi- 
ple phase slips case revealed the importance of the ratio 
u between the characteristic relaxation times of the am- 
plitude and the phase of the order parameter. While 
u « 1 is favoring simultaneous phase slips, evolutions 
with u » 1 more often happen with consecutive phase 
slips. Therefore, u » 1 is often a necessary condition 
for the relaxation to the ground state. The effect of the 
Langevin noise present in the equations is also studied. 
It appears that for higher values of the noise, the behav- 
ior may be different; but the general conclusions of the 
study are still valid. In some cases, phase slips may be 
favored by local inhomogeneities which could favor the 
simultaneous phase slips case, but this goes beyond the 
scope of the present study. 
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FIG. 5: Evolution of the distribution of the order parameter amplitude p (a)-(e) exhibiting the competition between simulta- 
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FIG. 6: Evolution of the distribution of the order parameter amplitude p in the presence of a Langevin noise of the order of 
1CT 1 . The phase slip are accelerated and happen less simultaneously. For it = 10 (d), we observe a mixing between consecutive 
and simultaneous phase slips. 



